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Overview 


This progress report for grant, NAG- 1-1 265, is given in the form of a manuscript which is 
currently under review for publication in the International Journal For Numerical Methods In 
Fluids, titled Sensitivity Analysis, Approximate Analysis, and Design Optimization For Internal 
and External Viscous Flows,” by Arthur C. Taylor, III, Gene W. Hou, and Vamshi Mohan Korivi. 
The manuscript was also presented as unpublished AIAA paper 91-3083 at the AIAA Aircraft 
Design Conference in Baltimore, MD, in September, 1991. The work illustrates the successful 
completion of all the tasks which were promised in the first year of the grant. 
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ABSTRACT 


Fundamental equations of aerodynamic sensitivity analysis and approximate analysis for the 
2D thin-layer Navier-Stokes equations are reviewed, and special boundary condition consider- 
ations necessary to apply these equations to isolated lifting airfoils on “C” and “O” meshes 
are discussed in detail. An efficient strategy which is based on the finite element method and 
an elastic membrane representation of the computational domain is successfully tested, which 
circumvents the costly “brute force” method of obtaining grid sensitivity derivatives, and is also 
useful in mesh regeneration. The issue of turbulence modeling is addressed in a preliminary 
study. Aerodynamic shape sensitivity derivatives are efficiently calculated, and their accuracy 
is validated on two viscous test problems, including: 1) internal flow through a double-throat 
nozzle, and 2) external flow over a NACA 4-digit airfoil. An automated aerodynamic design 
optimization strategy is outlined which includes the use of a design optimization program, an 
aerodynamic flow analysis code, an aerodynamic sensitivity and approximate analysis code, and 
a mesh regeneration and grid sensitivity analysis code. Application of the optimization method- 
ology to the two test problems in each case resulted in a new design having a significantly 
improved performance in the aerodynamic response of interest. 
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1.0 Introduction 


The tocus ot the present study is the continued development of efficient techniques for 
computing aerodynamic sensitivity derivatives for steady, viscous internal and external flows 
(Ret. [ 1 ]). In particular, this present work involves a direct extension of the recently developed 
methods ot Rets. [2] through [7], which have been successfully demonstrated on inviscid and 
viscous internal channel flows, to the classic external flow problem of an isolated airfoil. This 
extension to external flows is in essence a boundary condition problem. 

Sensitivity derivatives are defined as the derivatives of the system responses of interest 
(e.g., the specific thrust (F s ) of a nozzle, or the lift (C L ), drag (C D ), and pitching moment 
(Cm) coefficients of an airfoil) taken with respect to the design variables of interest (e.g., the 
parameters which define the geometric shape of the system, such as the thickness and camber of 
an airtoil). Once computed, these sensitivity derivatives are potentially useful in many ways. For 
example, the sensitivity derivatives can be used in approximate analysis, where if the changes 
in the design variables are small, resulting changes in a system’s response(s) can be accurately 
estimated, resulting in significant savings in computational costs. In addition, one of the most 
important applications of sensitivity derivatives is in engineering design optimization. The use of 
sensitivity derivatives in aerodynamic design optimization is demonstrated in the test problems 
to be presented. 

The present study is therefore design oriented, with the ultimate goal of the work being 
the development of tools which can be used by design engineers together with modern CFD 
(Computational Fluid Dynamics) software in improving the aerodynamic performance of the 
systems to which these codes are applied. Recent research efforts in the subject of aerodynamic 
sensitivity analysis with applications to design optimization have been intensive, and Refs. [8] 
through [20] is a representative (but not exhaustive) list of closely related works by other 
researchers, given here to provide additional background material on the subject, and its status. 

The remainder of the work is organized as follows: After the introduction, the next section is 
a presentation of theory, a section which is further sub-divided into six sub-sections, including 1) 
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governing equations, 2) spatial discretization and implicit tormulation, 3) fundamental equations 
ot sensitivity analysis and approximate analysis, 4) boundary conditions, 5) grid sensitivity, and 
6) ancillary sensitivity relationships. In the next section, the computational results are 2 iven in 
application to two test problems, including an internal flow through a double-throat nozzle, and 
an external flow over a NACA 4-digit airfoil. The final major section is a summary of the work 
where conclusions are given. 

2.0 Presentation of Theory 


2.1 Governing Equations 


The governing equations in this study are the 2-D thin-layer Navier-Stokes equations, given 


as: 


where: 


1 dQ 


™ Ot dr, + dr, 


( 1 ) 


( 2 ) 


Q = [p< pu, pv, ^>e 0 ] T (3) 

R(Q) is called the residual, and is clearly equal to zero for a steady-state solution. Q is a vector 
ot conserved variables, p is density, u and v are velocity components in Cartesian coordinates, 
and eo is total energy (i.e., eo=e + M t v > where e is the specific internal energy of the fluid). 
The inviscid flux vectors, F(Q) and G(Q) are given by: 

F(Q) = %F(Q) + %G(Q) 

J J (4) 

G(Q) = ^F(Q) + ^G(Q) 


4 



A transformation to generalized (£, 7 ) coordinates from Cartesian (x,y) coordinates has been 
made in Eq. ( 1 ), where are metric terms, and J is the determinant of the Jacobian 

matrix of this transformation. The Cartesian flux vectors F(Q) and G(Q) are given by; 

F(Q) = [pu- pu~ + P. /Hiv. (/)e 0 + P)u] T 

r T (5) 

G(Q) = [/n\ /mv, pv 1 + P. (/>e 0 + P)v] 

The pressure. P, is evaluated using the ideal gas law: 


P = 


:* - 1) 


pen - p 


^)] 


( 6 ) 


and ' is the ratio of specific heats, taken to be 1.4. The thin-layer viscous terms in generalized 
coordinates are given by: 


Gv(Q) = 


JL 

Rer 


[§ v 1 ’ gv 2 • gv, ! §V4 


where: 


gv, = 0 

gv, = + a 3 v n 

g v , = a 3 u n + aov v 
gvv = T Ql ( u2 ), + ^(v 2 ), 




i+i4 

). 0 = (s 

j 3 j 

/ ■ l> 

1 y\ 

7x + ny 

3 j r 

04 = J 


(7) 


( 8 ) 


The molecular viscosity is given by p, Stokes’ hypothesis for the bulk viscosity (A = 
-2p/3) has been used, a is the speed of sound, Pr is the Prandtl number (taken to be 0.72), 
and ReL is the Reynolds number. Nondimensionalization of Eq. ( 1 ) is with respect to p x 
and L'oc, the freestream density and velocity, respectively. The physical coordinates (x.y) are 
nondimensionalized by a reference length, L, and the viscosity is nondimensionalized by p^, the 


5 



molecular viscosity ot the freestream. The nondimensional molecular viscosity can be computed 
using Sutherland’s law and a reference temperature. T x , the static temperature of the freestream. 
For additional simplicity, however, in the laminar flow calculations of this work, the molecular 
viscosity is taken to be constant, equal to that of the freestream. For turbulent flow calculations, 
the algebraic turbulence model of Baldwin and Lomax [21] is used. 


2.2 Spatial Discretization and Implicit Formulation 


The governing equations are solved in their alternative integral conservation law form using 
an upwind cell-centered finite volume formulation. Only an overview of this method is presented 
here, with additional details found in Refs. [22] through [28], With this approach, the residual 
at each cell becomes a balance of inviscid and viscous fluxes across cell interfaces. As an 
example, this flux balance for the jk 115 cell in a typical computational grid is given by Eq. (9), 
for a steady-state solution, and for A£ = At; = 1 


Rjk - Fj + i - F;_i + G^j. - (X_i - G 


J-2 


'k +1 




+ G! 




(9) 


where subscripts j,k in Eq. (9) refer to the £,r? directions, respectively, and subscripts j ± ± 
refer to the £ = constant cell interfaces of the jk 111 cell, subscripts k ± i refer to the = 
constant cell interfaces ot the jk 1 * 1 cell. (All references to quantities which are evaluated at 
the cell interfaces will therefore require only a single subscript, either j ± \ or k ± I ) R jk is the 
discrete representation of the residual at the jk ^ 1 cell. Upwind evaluation of the inviscid fluxes is 
accomplished by upwind interpolation of the field variables, Q, from the approximate cell centers 
to the cell interfaces, where the flux vector splitting procedure of van Leer [29] is employed. A 
third-order accurate upwind biased inviscid flux balance is used in the stream wise (£) and in the 


normal ( 77 ) directions. The finite volume equivalent ot second-order accurate central differences 
is used tor the viscous terms. The resulting discrete higher-order accurate algebraic approximate 
representation ot the residual at each cell depends locally on cell-centered values of the vector 
Q at nine cells. That is, for the jk th interior cell 

Rjk(Q) = Rjk(Qj.k> Qj.k — 1 •» Qj.k+b Qj,k— 2» Qj,k+2 ’ Qj— l,k' Qj+l.k' Qj— >.k< Qj+2,k) 

GO) 
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Clearly adjustment are needed to Eq. (10) for interior cell equattons which are adjacent to 
the boundanes. When written for each cell (including boundary condition relationships to be 
discussed) and assembled globally, this can be expressed as 


{ R(Q* ) } = {0} 


where {Q } is called the “root” (i.e., the steady-state value of the field variables). Therefore. 
Eq. (11) represents a large coupled system of nonlinear algebraic equations, and thus finding 
a steady-state solution to the governing equations has been replaced (approximately) by the 
problem of finding the root, {Q*}, of this set of algebraic equations. 

The governing equations are discretized in time using the Euler implicit method, followed 
by a Taylor s senes hneanzation of the discrete equations about the known time level. This 
results in a large system of linear equations at each time step, given as 



1 1 


r<9R n (Q)i 

1 ) 

V 

LjAtJ 


!■ OQ j 

\) 


{"AQ} = R n (Q) 


( 12 ) 


{Q n+1 } = {Q n } + { n AQ} 
n = 1,2,3,... 


(13) 


Equations (12) and (13) represent the fundamental implicit formulation for integrating the 
governing equations in time to steady-state. In these equations, n is the time iteration index, and 
{ AQ} is the incremental change in the field variables between the known (n^) and the next 
(n+T h ) time levels. The matrix, [j^J, is diagonal, and contains the time term. 

In constructing exactly the true global Jacobian matrix, [- R q Q> ], of Eq. (12), both the 
interior cell equations as well as the boundary conditions must be considered, (although the 
contributions to this matrix from the boundary condition equations are typically neglected in 
many CFD codes). Considering the contribution to this Jacobian matrix from interior cell 
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equations only, the tour component vector equation which is associated with the jk th interior 
cell is isolated and extracted from the global linear system. Eq. (12). and is written below as 

[Aj{ n AQ J . k _ 1 } + [Bj{ n AQj k } + [C]{ n AQj k+1 } + [D]{ n AQj. k __>} + [E]{ n AQj, k+ ,} + 
[Fj{ n AQ J _,. k } + [G]{ n AQ j+1 . k } + [H]{ n AQj_, k } + [I]PAQ j+ _>. k } = { R]^ ( Q ) } (i 4) 

where { RjktQ) } is given by Eq. (10). and of course Eq. (14) represents the linearized form 
ot Eq. (10). The nine-point "difference stencil” represented by Eq. (14) is illustrated in Fia. 
(1), tor a typical interior cell. The nine 4x4 coefficient matrices [A] through [I] of Eq. (14) 

are constructed ot linear combination ot the inviscid and viscous flux Jacobian matrices, and 
[Bj also includes the time term. 

As a consequence ot Eq. (14), following global assembly of all interior cell equations, the 
resulting global coefficient matrix of Eq. (12) is sparse and has a banded structure, with nine non- 
zero diagonals, the individual elements of which are 4x4 block matrices. This matrix structure is 
illustrated in Fig. (2), which was taken from Ref. [30], Note that consistent implicit treatment 
ot the boundary condition equations and inclusion of these terms in the global coefficient matrix 
will sometimes severely disrupt the matrix structure which is illustrated in Fig. (2), depending of 
course on the type of boundary conditions. In addition to its use in Eq. (12) for time integration 
of the governing equations, this important Jacobian matrix, [^i], plays another central role 
in this study, which will be shown later. 

In principle, Eq. (12) can be repeatedly solved directly (using Eq. (13) to update the field 
variables), as the solution is advanced in time to steady-state, and for very large time steps, the 
direct method represents Newton's root finding procedure for nonlinear equations [31,32,33], 
The direct method however is not necessarily the most efficient procedure with respect to overall 
CPU time, and the large storage requirements of the method make its use not feasible in 3D. 
Therefore, more commonly, an iterative algorithm is selected for use in the repeated solution 
of Eq. (12). Popular choices of these iterative algorithms include approximate factorization 
(AF) [34], conventional relaxation algorithms [26,27,35], the strongly implicit procedure (SIP) 
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[36], and the preconditioned conjugate gradient method [37.38], to name a few. In the present 
research, the AF algorithm is used in the test problems to obtain steady-state numerical solutions 
to the governing equations. 


2.3 Fundamental Equations For Approximate Analysis and Sensitivity Analysis 

In this section, tundamental equations of aerodynamic sensitivity analysis and approximate 
analysis are reviewed, with additional details given in Refs. [1] through [7], Consider the 
vector, 3. the elements of which are independent variables which are typically called the design 
variables. Some, none, or all of the variables may be related to the geometric shape of the flow 
problem of interest, although the emphasis of the present study will be that of geometric shape 
variation. Computationally, the geometric shape of the domain is defined by the mesh upon 
which calculations are made, and the complete vector of (x,y) coordinates which defines the 
mesh is represented here symbolically as {X}. For a steady-state solution, the discrete residual 
vector given by Eq. (11) is rewritten in the following form 


{R(Q*(^),X(^),/i)} = {0} 


(15) 


where in the above, the direct dependence of the residual on the computational mesh, {X}, as 
well as its direct dependence (if any) on the vector 3 is now emphasized explicidy. Direct 
differentiation of Eq. (15) with respect to 3 k , the element of 3, yields 


dRl ( 3Q * j 
.^qJ l 03 k J 

Term 1 


'dRl f OX ] [OR] 

dxJ\d f y + Ud k / 

Term 2 Term 3 


(16) 


Equation (16) is an exact derivative of the discrete algebraic residual vector, and represents 
the central and most general relationship upon which those which follow in this section are 


based. The Jacobian matrix, 


dR 


, of Term 1 of Eq. (16) is identical to that found in the 


fundamental implicit formulation for time integration (Eq. (12)), and is evaluated here at the 
steady-state solution. It is thus well understood. The solution vector, ( ^ j, is the sensitivity 


of the complete vector of field variables with respect to the k 111 design variable. The matrix. 
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'or 

ox 


, ot Term 2 is the Jacobian ot the steady-state discrete residual vector with derivatives taken 


with respect to the complete vector of (x,y) grid coordinates, and is documented in Refs, [3,5]. 
The vector. {|J}, of Term 2 contains what is known here as the grid sensitivity terms, and 
is discussed in greater detail in Ref. [4], and also will be given special consideration later in 
the present study. The vector, [||-|, (Term 3) accounts for derivatives resulting from direct 
dependencies (if any) of the residual vector on 3 k . 


In the event that J k is not a geometric shape related design parameter (e.g.. the back pressure 
in a subsonic nozzle, or the angle of attack (a) or freestream Mach number (M-*.) for a airfoil) 
then Term 2 ot Eq. ( 16) will be zero. If 3 k is a geometric shape design parameter, its entire effect 
on the residual will typically be felt through the grid, and Term 3 of Eq. (16) will generally 
be zero. Therefore, Eq. (16) becomes 


_ ^ 

dQ 

when variation of geometric shape is not involved, and becomes 


dQ * ) = (OR) 
d,d k J I #4 J 


(17) 


'dRl fdQM 
.dQj \ d3 k } 



(18) 


when is a geometric shape design parameter. Eq. (18) represents the central relationship 
which is successfully demonstrated in Ref. [4] for accurately computing aerodynamic sensitivity 
derivatives with respect to variation of geometric shape. 

An approximate version of Eq. (18) which is useful in approximate analysis for estimating 
the steady-state solution changes which occur in response to small but finite geometric shape 
variations is given by 


'dR' 

.dQ, 


{AQ*} 



(19) 


where the approach represented by Eq. (19) is studied in detail in Refs. [3] and [5] for internal 
inviscid and viscous flows, respectively. In developing approximate analysis methods, other 
approaches of interest which might be used are 


79R' 

dQ| 


{AQ*} % 



( 20 ) 
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for variations other than 


geometric shape, (which is an approximate version of Eq. (17;), and 


'0R' 

. () Q. 


{AQ*} 


v;r 

dx 


ox 

d\ 


AJ k 


( 21 ) 


lor geometric shape variations only (which is a minor variation of Eq. (19)). 

In a typical design sensitivity analysis, the solution vector. j, of the preceding equations 
will provide far more information than is actually sought, and only a relatively small subset of 
this vector is extracted for use. For example, if the sensitivity derivatives of the aerodynamic 
torces on a solid surface boundary are to be calculated for a viscous flow, then the sensitivity 
derivatives ot the surface pressures and velocity gradients at the wall will be needed, which can 
be obtained as a subset of This is explained more fully in a subsequent sub-section, 

where the ancillary sensitivity equations of specific use herein are presented. 


2.4 Boundary Condi tions 

In the implementation of the fundamental equations of aerodynamic sensitivity analysis and 
approximate analysis, which were reviewed in the previous sub-section (Eqs. (15) through (21)), 
the consistent treatment and inclusion of the boundary conditions in these equations is essential, 
and must not be considered optional, as it typically is in the integration of the equations in 
time (Eq. (12)). The severely erroneous results which can arise as a consequence of failure to 
properly treat the boundary conditions are illustrated in an example problem in Ref. [2], with 
additional documentation on boundary condition treatment given in Ref. [6], 

For the isolated lifting airfoil problem, where the governing equations are typically solved 
on “C” or “O” meshes, the consistent treatment of the boundary conditions (during sensitivity 
analysis) is considerably more involved than that which is typically encountered in handling 
the boundary conditions for standard internal flow problems. Therefore, the objective of this 
section is to describe these difficulties, and to discuss various remedies. Many of the ideas to 
be presented in this section are taken from Refs. [33,39] involving the use of Newton (direct) 
solvers in obtaining inviscid and viscous steady-state solutions to airfoil problems. 
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2.4.1 Airfoil Surface Boundary Conditions 


The consistent treatment of the boundary conditions which are applied on the surface of the 
airfoil in this study (i.e., standard no-slip conditions) does not present any additional difficulty 
beyond that which is encountered for the boundary conditions which are used in typical internal 
flow problems. That is. the explicit application of the airfoil surface boundary conditions (as 
well as the application ot all boundary conditions which are used in the internal flow problem to 
be presented) can be represented at each point where the boundary conditions are applied (i.e.. 

at the boundary cell faces in the present study) as the solution of a boundary condition residual 
equation of the form given by 


{Rb(Qb(^)<Qi ( i) } = {0} 


( 22 ) 


where {R B } is a nonlinear tour component vector function of (at most) the local field variables, 

Q b* at the b0Undary celi face - Qf ’ the local ^ld variables at the first interior cell adjacent to the 
boundary, X,. the local coordinates of the grid, and also the very likely possibility of an explicit 
dependence on 3 is included. Differentiation with respect to 3 k , the k th design variable, yields 


f ^R b 




L<9QbJ \ 03], 
when geometric shape is not involved, and 


<9Rb 


3Qi 


l/£on _ 

Jl^kJ 


/3Rb 1 
l d3 k j 


(23) 


fdR 


B 


^QbJ 1 03 ], 


OQ 


OR 


B 


dQi 


fgQM _ dR B ] ( OX, 1 
1 03 k J OX i J 1 03], J 


(24) 


when geometric shape variables are involved. For approximate analysis, Eq. (23) can be written 
as 


<9R b 

L^QbJ 


{^Qb}- 


3R 


B 


l<9Qi 




if geometric shape is not involved, and Eq. (24) can be written as either of the following 
equations 


dR B 

L^Qb. 


{^Qb}- 


#r b 

aQi 


{^Qf} 


dR B 


OX, 




0 Rb 

<9Qb 


{^Qb) - 


0 Rf 


L3Qi 




0 R b 


OX, 


(-1 

1 d3 k j 


(25) 
two 

(26) 
(27) 
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for predictions of geometric shape change. 

Equations (23), (24), (25), (26), and (27) are of course fully a part of the global linear 
systems of equations given by Eqs. (17), (18), (20), (19) and (21), respectively, and can be 
explicitly included and solved simultaneously with the interior cell equations during the solution 
of these linear systems. Alternatively, yet equivalently, these boundary conditions equations can 
be pre -eliminated (i.e., pre-solved and substituted into the remaining interior cell equations) prior 
to simultaneous solution of the global linear systems. Pre-elimination of the boundary condition 
equations is the approach selected in the present study, and is explained in greater detail in Ref. 
[2], Note that in general, contributions to the Jacobian matrix, [§£], of the right-hand sides of 
Eqs. (16), (18), (19) and (21) can also occur (for some boundary condition types) from these 

boundary condition equations (in addition of course to the contributions to the matrix I— 1 ) 

’ L"QJ 

It is significant to note that simple boundary conditions of the form given by Eq. (22) result 
in consistently linearized boundary condition equations (i.e., Eqs. (23) through (27)) which, when 
included, do not in any way alter the basic structure of the coefficient matrix [||] , illustrated in 
Fig. (2) (i.e., the nine -diagonal, 4x4 block structure is preserved exactly), and thus direct solution 
of the fundamental systems of equations of sensitivity analysis and approximate analysis is not 
further complicated by inclusion of these boundary conditions equations. Unfortunately however, 
this is not the case in general, for all boundary condition types, as will be seen subsequently. 

2.4.2 Periodic Boundary Conditions For “C” and “O” Meshes 

When calculating flows over airfoils, in the likely event that a “C” or “O” type mesh has 
been selected for the calculations, '‘periodic” boundary conditions are applied along the “wake 
cut of the computational grid. There are different ways that these periodic boundary conditions 
could be applied. In the present research, the explicit application of these boundary conditions 
is not expressed (as before) as the solution of special boundary condition equations of the type 
given by Eq. (22). Instead, interior cell residual equations for the interior cells immediately 
adjacent to the periodic boundaries (i.e., the wake-cut) are each expressed as functions of 
cell-centered values of the vector Q at nine cells in the domain (i.e., for a higher-order accurate 
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upwind spatial discretization), where this functional relationship in physical space is identical 
to that expressed by Eq. (10) tor a general interior cell. This means that wherever necessary, 
the evaluation ot the inviscid flux vectors ot Eq. (9) is accomplished by a consistent upwind 
interpolation of the field variables from local interior cells across the wake cut. and that the terms 
ot the viscous flux vectors are evaluated using consistent central differences across the wake cut. 

Although the tuncfional relationship expressed by Eq. ( 10) for a general interior cell equation 
is preserved in physical space tor interior cells involving periodic boundary conditions, clearly 
it is not preserved for these cells (with respect to the structure and ordering of the cells) in 
the computational domain. That is, when periodicity is involved, the interior cell equations 
depend explicitly on the field variables at local neighboring cells in the computational grid, and 
in addition, become functionally dependent on the field variables of cells which are quite distant 
in the computational ordering. Consequently, the linearized version of Eq. ( 10), given by Eq. 
(14) for a general interior cell, must be modified for interior cells involving periodicity, in order 
to properly account tor these boundary conditions. Periodic boundary conditions affect each first 
and second interior cell equation immediately adjacent to the boundary where “periodicity" is 
enforced, resulting (for the higher-order spatial discretization) in two periodic 4x4 matrix terms 
contributing to the left-hand side ot Eq. (14) for the first interior cell equation, and one such 
periodic term for the second interior cell equation. 

The most significant impact of the periodic boundary conditions is that there is some 
restructuring ot the coefficient matrix, |^j , where the neat, banded, nine-diagonal structure 
which is illustrated in Fig. (2) is no longer preserved exactly. This is illustrated in Ref. [33], 
and is ot course is a consequence of the previously discussed adjustments to Eq. (14) which 
are required for the interior cell equations involving periodicity. Furthermore, depending on the 
direction which is selected when ordering sequentially the individual cell equations for assembly 
into the global linear system (either the tangential (J) to the airfoil direction, or the normal 
(K) to the airfoil direction), the non-zero contributions to the global coefficient matrix from 
the periodic boundary conditions will fall either inside or outside of the central bandwidth of 
the matrix. (Note: The central bandwidth here refers to all of the matrix elements, either zero 
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or non-zero, which fall between the outermost diagonals, [H] and [I], of Fig, (2) and is thus 
the bandwidth ot the matrix, neglecting periodic terms, if any, which fall outside.) For “C” or 
"O mesh calculations, it the equations are numbered in the tangential (J) direction, the periodic 
coefficients will all fall inside the central bandwidth, but will all fall outside the central bandwidth 
(some at extreme distances) if the numbering is in the normal (K) direction. (Exception: There 
exists a special method explained in Ref. [33] for numbering the equations in the ”K" direction 
lor "C ' meshes where the periodic coefficients all fall completely inside the central bandwidth. 
This method results in a doubling of the central bandwidth of that which is achieved with a 

standard K ordering ot the equations, and hence this non-standard ordering was rejected in 
the present study). 

A J direction ordering ot the equations will therefore in principle allow the use of a pure 
banded matrix direct solver solution procedure which takes advantage of the fact (in terms of 
storage and work) that outside of the central bandwidth, all of the matrix elements are zero. In 
contrast, for a pure direct solution of the linear system, a standard K” direction ordering will 
require the use ot a lull matrix solver to account for the periodic terms, which is not feasible 
tor practical fluids problems involving the tull governing equations, even in 2D. 

For typical airfoil calculations, the “J” dimension of the grid is usually significantly larger 
than the K dimension. A standard “K” ordering of the equations will thus result in a 
significantly smaller central bandwidth ot the coefficient matrix than will a “J” ordering. To 
overcome this dilemma, a hybrid direct solver / conventional Richardson type relaxation strategy 

is proposed and implemented in the present study, for airfoil problems on “C” and “O” meshes, 
as follows: 


1) A K direction (i.e., normal to the airfoil) numbering of the equations is used in 


constructing the coefficient matrix, 
bandwidth. 


dR 


, in order to minimize the width of the central 


2) The coefficient matrix is split into two parts, such that 

<9R 


dQ 


= [M] + [n; 


(28) 
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where the matrix [M] contains all of the elements which fall inside of the central 
bandwidth of . Thus [M] has the matrix structure illustrated in Fig. (2). The 
matrix [N] is thus very sparse, and contains all of the non-zero contributions to 

[a Q j 

which tall outside the central bandwidth, resulting from the periodic boundary conditions. 

3) The matrix [M] is LU factored directly using a conventional banded matrix solver, to 
yield 

[M] = Mn (29) 

4) A conventional Richardson type iterative strategy is invoked which in principle could 
be applied to Eqs. (16). (17), (18), (19), (20), and/or (21). For example, if Eq. (18) is 
selected, the iterative strategy is 



i = 1 , 2, 3. ■ • ( imax) k ^0) 

k = 1, ndv 

where in the above, i is the innermost index, (imax) k is the maximum number of iterations 

required to converge the k th linear system to the desired tolerance, V is the outermost index, 

• . f 1 (* max ) k 

which associates the solution vector, with a particular k th design variable, and 

ndv is the total number of design variables for which sensitivity derivatives are required. It is 
emphasized that since [L] and [U] (in addition to and [N]) are constant with respect to 
the indices T and ‘k\ they only need be computed once and repeatedly reused. Following the 
LU factorization, a single iteration of Eq. (30) is inexpensive, of course, since [L] and [U] are 
lower and upper triangular matrices, respectively. 

The coefficient matrix, fq- , is diagonally dominant for a first order upwind spatial 
discretization, but is not diagonally dominant for a higher-order spatial discretization [26], and 
therefore, convergence of the iterative strategy represented by Eqs. (28), (29) and (30) is not 
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assured [35]. In the airfoil test problem to be presented, the iterative strategy of Eq. (30) was 
at first divergent, but was subsequently made convergent through the simple use of a single 
conventional under-relaxation parameter. (omega). Other remedies which are suggested as 
possible methods to correct the failure of this iterative strategy to converge include, but are not 
limited to: 

1) Use ot a first order accurate upwind spatial discretization for the inviscid terms of the 
first row of interior cells immediately adjacent to the wake cut. which will result in 
only a single non- zero periodic 4x4 matrix coefficient contribution (instead of three such 

coefficients) to the matrix [N] of Eqs. (28) and (30), for each point where periodicity 
is enforced. 

2) Use of a “ghost” cell equation at each point where periodicity is enforced, which is 
retained explicitly in the global coefficient matrix, resulting in a more strongly implicit 
treatment ot the periodic terms, with complete preservation of higher-order accuracy 
across the wake cut. This is the approach of Ref. [33], where difficulties with 
convergence of the iterative strategy were not encountered [39], 

3) Use of matrix pre-conditioning. 

In addition to the contributions to the matrix, , resulting from the application of periodic 
boundary conditions, there are also contributions to the matrix, |£- , of Eqs. (16), (18). (19), 
and (21), which must be considered tor viscous flow. For inviscid flow, there are no additional 
adjustments needed to to account for periodicity. However, a consistent treatment of all 
viscous metric terms in the application of periodic boundary conditions will result in the need 
for some adjustments to some of the terms of [§|] to properly account for the periodicity. A 
close examination of the documentation given in Ref. [5] on the details of the construction of 
the viscous terms ot will provide additional explanation and verification of this. 
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2.4.3 Far-Field Boundary Conditions 

The far- fie Id boundary conditions which are used in this study make use of the Riemann 
invariants, and these boundary conditions are given in Refs. [33] and [40]. For lifting airfoils, it is 
well known that significantly improved computational accuracy can be achieved with the addition 
of a "point-vortex” correction to these boundary conditions. This "point-vortex" correction is 
developed and presented in detail in Ref. [40], and will be addressed in the context of its treatment 
in sensitivity analysis in this sub-section. It this point-vortex correction is not included, then the 
explicit application of the tar-field boundary conditions at each location can be expressed exactly 
as the solution of a boundary condition residual expression of the form given previously by Eq. 
(22). In this event, the contribution trom the far-field boundary conditions to the equations of 
sensitivity analysis and approximate analysis (i.e., to Eqs. (16) through (21)) is therefore very 
straightforward, using Eqs. (23) through (27). 

The use of the "point- vortex” correction described in Ref. [40] to improve the far- 
field boundary conditions is straightforward to implement in an explicit sense. Its explicit 
implementation involves the use of a point-vortex (centered at the quarter chord) representation 
of the airfoil, where the strength of the point-vortex (i.e.. the circulation, D is proportional to 
the lift coefficient, C L , of the airfoil. The purpose of this point-vortex is to more accurately 
model the influence of the lifting airfoil on the velocity field in the vicinity of the far-field 
boundaries, (compared to the alternative of assuming a freestream velocity field here), resulting 
in more accurate airfoil calculations, particularly as the extent of the far-field boundary from 
the airfoil is decreased. 

The implementation of this point-vortex correction results in a numerical coupling between 
the far-field boundary condition equations and (through the lift coefficient, C L ) the field variables 
(and also the x,y grid coordinates) on and immediately adjacent to the surface of the airfoil. As 
a consequence of this coupling, there are algebraically messy, complex additions to the global 
Jacobian matrix, (which destroy the banded matrix structure, illustrated in (Fig. (2)) and 
also to [-jxj- In order to avoid the task of explicitly deriving these terms and their precise 
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locations in these global Jacobian matrices, a simplifying strategy is proposed, which makes use 
of the following ideas: 

1. In the discrete system of algebraic equations which models the steady-state fluid flow 
(i.e., Eq. (11)). the lift coefficient. Cl, is to be treated as an additional field variable. 

2. The chain rule is employed judiciously. 

3. In a typical sensitivity analysis for airfoil design, the sensitivity of the lift coefficient 
with respect to the design variables will be sought. Explicit expressions are therefore 
available for use in calculating the sensitivity derivatives of C L with respect to the design 
variables, following solution of the global sensitivity equations. (This expression for C L 
is given in a subsequent sub-section). 

4. Recall that a strategy which involves iteration will be used to solve the equations of 
sensitivity and approximate analysis to account tor the periodic boundary conditions 
(regardless of whether the point-vortex correction is used or not). 

When the point-vortex correction is included, the boundary condition residual equation which is 
solved at each far-field boundary cell tace can be expressed in a form similar to that of Eq. (22), 
except that an explicit dependence ot the boundary condition equation on Cl is now identified. 
That is, Eq. (22) becomes: 


{RB(Q|(j).Qf(,i).x,(j). j.c L )} = {o) (3i, 


where all the terms of Eq. (31) are as previously defined. Differentiation with respect to J k , 
and applying the chain rule on the term involving Cl, the result is 


<9Rb 

<9Qb. 



[0Rb| fdQM 

. dQl J t ddy } 


'dRel/dX,! fdR B ] f (9R b j dC L 
.dX,JU4J l<94 J + Uc: L /d4 


(32) 


The four component vector, j j , is very straightforward to derive analytically, since the 
explicit dependence of Rb on Cl is known and uncomplicated [40], There are also additions 


to the matrix. 



, which are straightforward to derive, which arise as a consequence of this 
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boundary condition correction. The scalar term, is evaluated using an expression of the 
form 


cICl 

c.U k 


r)C\ 


7£Q1 

/ \ <94 


+ 


*4 

dx 


ox ) or . 


+ 


r94 j <94 


(33) 


where the global column vectors, j ^ j and are sparse. Since the value of ^ is 

likely to be at least one of the desired final results of an airfoil sensitivity analysis, the precise 
terms involved in Eq. (33) (i.e., {^} T and {^} T ) have been derived analytically for this 
purpose, and are given in a subsequent sub-section. Note in Eq. (33) that the notation for a total 
derivative has been used on the left-hand side of the equation, indicating that the total rate of 
change is included in the expression, and thus distinguishing it from the partial derivative (which 
accounts for explicit dependencies of C L on 3 k , if any) on the right-hand side of the equation. 


It should be understood, however, that ^ is still a partial derivative in the sense that C L is in 
general a function ot multiple independent design variables. 

Specialization ot Eq. (32) to the case of either geometric or non-geometric shape design 
variables and to the case ot either sensitivity analysis or approximate analysis is analogous to 
that shown in Eqs. (23) through (27). In particular, for approximate analysis, the term ^ of 
Eq. (32) becomes AC L , and Eq. (33) becomes: 


AC l % 


d£i 

dQ 


T 

{^Q*} + 


dC L 

dX 


T 



(34) 


At this point, the tar-field boundary condition equations given by Eqs. (32) and (33) are in a 
convenient form for assembly into the global system of equations for subsequent solution in 
the hybrid direct solver / iterative strategy previously described, where iteration is also used to 
account for the periodic boundary condition terms (whether or not the point-vortex correction 
is included). Global assembly of Eq. (32) into the iterative strategy of Eq. (30) results in the 
following modified iteration strategy to account for the point-vortex correction 



i = 1, 2, 3, • • •. ( imax) k (35) 


k = 1, ndv 
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where in Eq. (35) the only non-zero contributions to the global vector. are trom the 

tar-field point-vortex corrected boundary condition equations. The scalar term is repeatedly 
calculated tor all 'i' and 'k' indices using the expression 


dC L 

<14 


0£i\ W) 

(9Q j \ c)J k J 


1-1 


+ 


d£i 

ox 


ox 


+ 


OCX 


rJJ k J 03 k 


(36) 


The iterative strategy of Eq. (35) thus avoids the need for determining explicitly the algebraically 
messy contributions of the point-vortex corrected boundary conditions to the Jacobian matrices. 

and [ox/ 


2.5 Mesh Sensitivity Analysis 


In this section, calculation of the terms of the mesh sensitivity vector, j j, of Eqs. (16), 
(18), and (21) is given special consideration. One approach which can be used for calculation 
ot these terms is the ‘brute force finite difference method, where the geometric shape design 
parameters are perturbed slightly, one at a time, and the mesh generation program (which was 
used to generate the initial mesh) is repeatedly re-run to generate perturbed grids. Then, each 
grid sensitivity vector is calculated as 



where in the above, central finite difference approximations are used for greater accuracy, but 
require the generation of two perturbed meshes per design variable instead of the one which 
is required using a forward finite difference. The principle advantage of this approach is that 
it is conceptually simple and easily applied. For this reason, it is also well suited for use 
in an automated design optimization loop. Since by design the equations of mesh generation 
are typically very smooth, the method can be expected to be reliable in producing satisfactory 
accuracy. Furthermore, tor simple 2D geometries using algebraic grid generation, the procedure 
should be computationally efficient, even when the number of design variables is large. 

High quality elliptic and hyperbolic grid generation codes, which are often selected for 
generating “C and "O” airfoil meshes, are far more computationally expensive to use than are 
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their algebraic counterparts, and thus use of the “brute force” method to obtain grid sensitivities 
in this case could become excessively costly, particularly if the number of design variables is 
large, and it the computations are to be made repeatedly in an automated design loop. One 
strategy proposed in Ret. [4] to help restore computational efficiency to this problem would be 
to differentiate analytically the equations of mesh generation (for the particular mesh generation 
code of choice), to determine the Jacobian matrix of the entire grid, X, with respect to the grid 

points on the boundary of the domain, X s . Then grid sensitivity derivatives could be efficiently 
calculated using the relationship 



Other researchers [41] have implemented the approach of Eq. (38) in application to the airfoil 
grid generation program of Ref. [42], Another important consideration with respect to some 
sophisticated grid generation codes is that they are used interactively, which would prohibit their 
use (interactively) in an automated design loop. 

In order to help overcome these difficulties which have been discussed, a procedure is 
proposed herein for use in the repeated evaluation of grid sensitivity terms, as well as for use 
in efficient grid regeneration. The method, which will be presented subsequently, is based 
on an “elastic membrane” analogy to represent the computational domain, and grid sensitivity 
derivatives are calculated from a standard structural analysis computer program using the finite 
element method. 

As the geometric shape of the flow domain is continuously changed, as required by any 
geometric shape optimization process, the mesh points in the domain must be properly adjusted 
m the design iterations, to avoid the numerical errors induced by excessive mesh distortion. 
The requirement of mesh regridding distinguishes shape design optimization from other design 
optimization applications. A simple way for automatic mesh regridding can be established by 
introducing a set of basic displacement vectors, V k , to describe the patterns by which the mesh 
is regridded. The relationship between the original mesh, X Q , and the regridded one, X, can 
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then be expressed in the torm ot a linear combination of these basic displacement vectors, and 
their associated weighting coefficients, J k , as 

ndv 

X = Xq + y, J k V k (39) 

k = 1 

where in the above, ndv is the number of design variables, if the weighting coefficients are taken 
to be the design variables. The vector. X 0 , represents the initial mesh, which is produced using 
uny conventional mesh generation code ot choice. In this case, the basic displacement vector, 
V k , is simply equal to the required mesh sensitivity vector, j|^j, of Eqs. (16), (18). and (21). 
That is, the grid sensitivity derivatives are calculated by differentiation of Eq. (39). which yields 



Note that the grid sensitivity vectors, { V k }, do not change when the design variables are changed, 
provided that the domain is always regridded using Eq. (39), as the shape of the domain changes. 
Therefore, these grid sensitivity derivatives need only be calculated once and then stored prior 
to the start of an aerodynamic optimization strategy, and they can be repeatedly re-used as often 
as needed tor grid sensitivity analysis, as well as for automatic mesh regeneration. 

The basic displacement vectors, V k , can be in any form; however, they must be independent 
ot each other. In structural shape design optimization, the elastic displacements induced by the 
boundary perturbations are commonly selected to represent the basic displacement vectors. In this 
way, the movement ot the mesh points is governed by the nature of linear elasticity, which not 
only preserves the continuity of the mesh, but also avoids any mesh-overlapping. It is proposed 
herein that the same practice can be applied to aerodynamic shape optimization problems, in 
which an imaginary elastic medium is introduced to represent the computational domain. 

More specifically, the basic displacement vectors can be generated by either the fictitious 
load method [43], or the prescribed displacement method [44], The former produces the basic 
displacement vectors by applying a unit load at each of the nodes along the boundary in the 
direction along which the node is allowed to move. This concept is illustrated in Fig. (3), for 
a representative airtoil grid. The latter, however, produces the basic displacement vectors by 
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imposing a non-zero displacement (in response to a unit change in each design variable) along 
the varied boundary. This concept is illustrated in Fig. (4), for a representative airfoil grid. 
The fictitious load method is usually applied to the case where the location of each node on the 
v aried boundary is considered as a design variable, whereas the prescribed displacement method 
is_good in the case where the shape of the boundary to be designed is parameterized 

In the following, a NACA 4-digit airfoil is employed as an example to demonstrate the 
application of the prescribed displacement method for mesh-regridding in an aerodynamic shape 
optimization environment. The profile of the NACA 4-digit airfoil can be precisely represented 

by polynomials in terms of the maximum thickness, T. the maximum camber. C, and the location 
of maximum camber, L, as 


( _ f f < x ) + C(2Lx - x 2 )/L 2 , x < L 

lf(x) + C(l-2L+2Lx-x a )/(l-L) a . x>L <4I) 

where 

t(x) = ±.5T(0.2969\/x - 0.126x - 0.3516x 2 + 0.2843 x 3 — O.lOlox 4 ) (42) 

where, the ± in the expression tor f(x) is V for the upper surface of the airfoil, and for 
the lower surface. 


Since the denvatives of the airfoil shape with respect to T, C, and L are continuous, it 
is understood that small changes in T, C, and L will induce small changes in airfoil shape. 
Theretore, according to a Taylor’s series expansion, such a change in the airfoil shape can be 
expanded approximately into a linear function of AT, AC, and AL, given as 


y(x) = y 0 (x) + 


dyo(x) 

dT 


AT + 


dyo(x) 

dC 


AC + 


dy Q (x) 

dL 


AL 


(43) 


where 

AT = T - T 0 

AC = C - C 0 (44) 

AL = L - L 0 

where T 0 , C 0 , and L 0 are the initial values of these three shape parameters associated with the 
initial airfoil shape, y 0 (x), and the initial grid, X 0 . 
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The derivatives ' - ^ x ■ , dy ^\ and 
the allowable changes in the airfoil’s shape, 
by Eq. (39) as 


in Eq. (43) represent special patterns which control 
The new mesh, X can be defined in a form given 


X = X 0 + AT • V, + AC • V, + AL • V 3 (45) 

where AT, AC, and AL are taken to be the design variables for equivalently T, C. and L 
are the design variables, through Eq. (44)). The basic displacement vectors, V,, V 2 . and V 3 
can be obtained by the prescribed displacement method, discussed above. These vectors are 
obtained numerically through implementation of a finite element model, with each cell in the 
computational mesh being considered as a plane stress quadrilateral element. A finite element 
matrix equation is formed to solve for each basic displacement vector (i.e., the movements of 
all the grid points) which is realized throughout the elastic membrane model of the domain, in 
response to the non-zero boundary movement which is specified through Eq. (43), for a unit 
change (or some other conveniently scaled change) in each of the design variables. Note that 
the finite element matrix equation is linear with a symmetric and banded coefficient matrix, 
which is solved directly by a single LU factorization, followed by multiple re-uses of this LU 
factorization as multiple solutions (i.e., one solution, V k , for each design variable) are obtained 
for multiple “load vectors” through simple efficient forward and backward substitutions. 

It is clear that Eq. (43) represents a particular parameterization of the surface of the 
airfoil which will only closely approximately the NACA 4-digit parameterization defined by 
Eqs. (41) and (42) if AT, AC, and AL are small. However, if during design it is not a 
requirement to remain exactly within or close to the allowable shapes defined by the NACA 
4-dign parameterization, then of course Eq. (43) is a valid (but different) parameterization of 

the airfoil shape, even for large AT, AC, and AL. Thus the classic NACA 4-digit airfoil is 
presented here only as an example. 

In order to demonstrate and validate the mesh regeneration strategy, a numerical example is 
given. Starting with an extremely coarse initial grid, X c , with only 21x10 points (for illustration 
purposes), for a NACA 1412 airfoil (i.e., T o =0.12, C o =0.01, L o =0.40), generated using the grid 
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generation code of Ref. [42], the domain was regridded using Eq. (45) for a new airfoil shape 
defined by Eq. (43) and (44), for C=0.04, T=T 0 , and L=L 0 . The initial airfoil and grid is shown 
in Fig. (5a). where in Fig. (5b), the new airfoil and grid are shown. The changes which are 
seen in the new mesh resulting from the application of this methodology appear satisfactory in 
this test example. 


2.6 Ancillary Sensitivity Relationships 


The purpose ot this sub-section is to include some important extra terms and relationships 
which are used in the present study in computing sensitivity derivatives, in order to make the 
presentation of the methods more complete. Specifically, expressions are given for generalized 
aerodynamic force and moment coefficients in 2D, and their sensitivity derivatives. In addition, 
expressions are derived for use in computing the sensitivity derivatives of the thrust, mass flow 
rate, and specific thrust of a nozzle. 

Fig. (6) illustrates the i th element (oriented at an arbitrary angle in space) which is located 
on the boundary of the geometric shape of interest, over / through which the fluid is passing. 
In the figure, the coordinates (x bi ,y b .) and (x bi+l ,y bi+1 ) are the physical (x,y) coordinates at 
either end of this i m element, and are assumed to be nondimensionalized by L, the reference 
length. The convention is established that as one moves along the surface in the direction of 
increasing the index, *i\ then the solid surface is on the right, and the fluid is on the left. 

Nondimensional pressure (C p ) and skin friction (Cf) coefficients which are associated with 
this i* element on the boundary are defined as follows 


C P , = 


O Poo U 5c 


Cf,= 




2PocX T ^o 


(46) 


where in the above, P bi and are the pressure and shear stress, respectively, which are 

associated with the i* element on the boundary, and is the dynamic pressure of the 
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treestream. Nondimensional force coefficients in the x and y directions, given as C x , and C’ v , 
respectively, for the surface element, are given as 

( x > ~ ( Pi (>'bi+i _ y b; ) + ( ’r, (x bi+1 - x bi ) 

,, r , , , (47) 

( -y. ~ '-pi ( x bi - x b i+ , ) + C fi (y bl+1 - v bi ) 

The total force coefficients in the x and y directions are of course given by summing the above 
expressions over the total number of elements of interest, NE, and the result is 


NE 

( x = ^ C x , 

NE (48) 

<y = £c„ 

i = l 


Sensitivity derivatives ot these force coefficients taken with respect to 3 ^ are given as 


dC’ x 

d/i k 


NE 



~ x bj ) + Cfj 


V d3 k 


dxbj 

d3 k 


dC’y 

d4 



d£p, 

ds k 


( x bi ~ x b i + l ) + 





y bi ) + Cfj 


dyb i+ , 

d3 k 



(49) 


Note in the above expressions that terms such as Jg-, etc., are evaluated as being elements 
ot the vector j^-J of the previously discussed right-hand side of Eq. (16), (18), or (21) and 
and also Wt ^ obtained from the vector { } by solution of these equations. (Note: 

Since Cf, involves gradients of the velocity at the i 1 * 1 element of the boundary, then 
involves terms from both j^-j and resulting in an expression not shown here which 

is algebraically complex.) 

If the sensitivity derivatives of the lift (C L ) and drag (C D ) coefficients of a body (e.g., an 
airfoil) are desired, it is only necessary to convert to a new coordinate system with axes aligned 
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tn the L and D directions, which are rotated at an angle of attack, a, with respect to the x and 
y axes. The litt and drag coefficients are then calculated as 


C L — f v cos o — ( ’ x sin a 


(50) 


Co = C v sin a + C x cos a 


(51) 


The corresponding expressions for the sensitivity derivatives of the lift and drag coefficients are 


<ic L 

dC v 

1 dC x 

d,i k 


cos q — - sin a 

d.4 

dC D 

dC’y 

dC x 

cU k 

~ d4 

Sin Q + — — - COS Q 
d^k 


(52) 


(53) 

Note Eqs. (64) and (65) are not valid if the angle of attack, a, is the design variable (i.e.. if a = 
4). In this special case, these expressions would simply have additional terms. (See Eq. (33).) 

In addition, for airfoil calculations, the sensitivity derivatives of the pitching moment 
coefficient, C M , could likely be of interest in a design problem. The moment coefficient. C Mi , 
which is associated with the i* element (see Fig. (6)) on the solid wall boundary of interest, 
about the point (x 0 , y 0 ), with the convention that a positive moment is clockwise, is given by 

C .Mj = C Xi (y Ci - y 0 ) - C Vl (x Ci - x 0 ) (54) 

where the coordinate (x Cj ,y Ci ) is the midpoint of the i* element on the boundary. That is 

X Q - 2 ( X bi +1 + x b,) y Ci = 2 (ybi +1 + ybi) (55) 

The total moment coefficient is of course the sum of all the elemental moment coefficients over 
the total number of elements on the boundary of interest, given by 


NE 

C’m = C.Wi 

i=i 


(56) 
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The sensitivity derivative of the pitching moment coefficient 
variable is thus 



with respect to the k ,h design 


x o - xj + Cy, 




(57) 


where 


C ^L = 1 ( r)x ' b ' + ' ■ ) foo _ 1 / r) y bj+1 r7y bi \ 

^ 2 v rTi k ; a j k - 2 ( m~ + OA ) 


(58) 


In the present study, an internal flow problem was considered where the sensitivity derivatives 
of the thrust, T, mass flow rate, m, and specific thrust. F s , of a nozzle are of interest, and the 
xpressions used in these calculations are developed next. The total thrust which is produced by 
a nozzle can be calculated from the previously defined force coefficient, C* , summed over all 
of the interior solid walls, plus an additional term involving the integrated static pressure over 
the inflow boundary of the nozzle. That is 


NEJ 

T = + zl C Pj ( - v j+! - - v j ) (59) 

j=i 

where m the above, a positive thrust acts in a direction opposite to the positive direction of the 
x axis, the y axis is assumed to be in or parallel to the plane of the inflow boundary, C' Pj is the 
static pressure coefficient of the j' h element on the inflow boundary, yj and y j+1 are the starting 
and ending coordinates, respectively, of the j' h element on the inflow, where the direction of 
increasing the index, ‘j\ and the positive y direction are in this case defined to be always the 
same, and NEJ is the total number of elements on the inflow boundary. Nondimensionalization 
of T in Eq. (59) above is of course the same as defined previously for C x . (Note: The expression 
above does not consider the aerodynamic forces, if any, due to flow over the exterior surfaces 
of the nozzle.) The sensitivity derivative of the thrust is 
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The mass flow rate, rh, through the nozzle can be computed using the expression 

NEJ 

m p i UjLvj+t -.Vj) ( 61 ) 

j=i 

and the sensitivity derivatives become 



Often with nozzle calculations, the specific thrust, F s , is a system response of interest, and is 
defined by 


F s = 


T 

rh 


(63) 


and the sensitivity derivatives can be calculated using expressions previously defined 

JF. 

d.i^ rh- 


(64) 


3.0 Computational Results 

3.1 Internal Flow — Double-Throat Nozzle Problem 

3.1.1 Description of the Test Problem 

The first test problem to be presented is that of an internal flow through a double-throat 
nozzle, where the flow is accelerated from a Mach number on the inflow boundary of about 
0.10, to a Mach number which exceeds 2.80 at some places on the outflow. The Reynolds 
number, RE L , is 100, based on a reference length, L, of one-half the height of the nozzle at 
the smaller of the two throats. Other researchers have conducted numerical studies on the flow 
through this nozzle geometry, with documentation provided in Refs. [5,45,46,47], 
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A computational grid is used with 171 points evenly spaced in the streamwise direction, 
and 38 points with grid stretching in the normal direction, to resolve viscous gradients in the 
vicinity ot solid walls. Fig. (7) illustrates the grid and geometry. Boundary conditions are 
pecified as follows: 

1) On the inflow, entropy and stagnation enthalpy are held constant at the freestream values, 
the v component of velocity is zero, and the u component is extrapolated from the 
interior of the domain. 

2) On the outflow boundary, all variables are extrapolated. 

3) On the lower wall, velocity no-slip is enforced, and the wall static temperature is specified 
to be that of the stagnation temperature of the freestream. 

4) On the upper boundary of the computational domain (i.e., along the centerline of the 
nozzle), flow symmetry boundary conditions are enforced. 

The Mach contours tor the steady-state laminar flow solution through the nozzle are shown in 
Fig. (8). 

The geometry of the nozzle is defined parametrically using analytical expressions which are 
presented subsequently, where this material is also given in Ref. [47], Because of symmetry, 
it is only necessary to parameterize either the lower or the upper solid wall surface. The upper 
wall is selected here. The wall is described using five polynomial arcs, as illustrated in Fig. 
(9). In all cases, continuity of slope is enforced at the transition point from one arc to the next. 
Continuity of curvature is also preserved at the transition point from arc II to arc III and from 
arc III to arc IV, but is discontinuous as the transition from arc I to arc II and from arc IV to 
arc V. Using Fig. (9), these five arcs are defined as follows 


ARC III (X 3 < X < X 4 ) 


Y = H + f 4 1 X 


X 2 

6 


X 3 + X 4 


X + X 3 X 4 


(65) 
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From Eq. (65), the following is deduced 

Y(Xi) = H + (f )<X:» , (2X 4 - 

Vx(X:, ) = (t) (X - i) 2 ( Xi_ t) 

/ A \ / Y \ (66) 

Y(X 4 ) = H + f-J(X 4 )M x 3 - y ) 

Y<(X 4 ) = |(X 4 ) 2 ^X 3 -y) 

where in the above and henceforth the subscript x indicates differentiation with respect to X 
(i.e., \ x = jy)- 


ARC II (X 2 < X < X 3 ) 


Y = Y(X 3 ) + Y x (X 3 )(X-X 3 ) 


1 - - 



X - x 3 
x 2 - X, 


(67) 


From Eq. (67), the following is noted 


Y(X 2 ) = Y(X 3 ) + ^(X 2 - X 3 )Y x (X 3 ) ( 68) 


ARC I (Xi < X < X 2 ) 
Y = Y(X 2 ) 


(69) 


ARC IV (X 4 < X < X 5 ) 


Y = Y(X 4 ) + (X - X 4 ) Y x (X 4 ) + B(Zi) 2 -h C(Zi; 


(70) 


where 

Z 1 = (X-X 4 )/(X 5 -X 4 ) 

B =4D — 3Y x (X 4 ) 

(71) 

C = —3D + 2Y x (X 4 ) 

D = [Y(X 5 ) - Y(X 4 )]/(X 5 - X 4 ) 

Note that X5 is the location of the second throat, and Y(Xs) represents the half-height of this 
second throat. 
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ARC V (X 5 < X < X 6 ) 

V = V(X g ) + [Y(X 6 ) - Y(X 3 )](Z 2 ) 3 (2 - Z>) 


( 72 ) 


where 

z_, = (X - X 5 )/( X 6 - X 5 ) (73) 

As a consequence of the above relationships, there are ten geometric shape parameters which can 
be identified from these expressions, which in the present work will also detine the elements of 
the vector 3 . known here as the design variables. These parameters are given below, including 
their numerical values for the “baseline” test geometry (i.e., the “initial design” of the nozzle). 


where 


i • 3-2 , 3i 

• 3^ , 3s , 3s i 

3j . 3 $ , 3$ 

■h = 

H 

+ 1.0 

3-2 = 

X 3 = 

-4.0 

3;i = 

x 4 = 

+2.3 

3 4 = 

A 

-0.03 

3s = 

X, 

- 10.0 

3e = 

Xi = 

- 12.0 

3i = 

X 5 = 

+7.0 

3n = 

Y(X 5 ) = 

+ 1.6 

139 = 

X 6 = 

+ 14.0 

.^10 = 

Y(X 6 ) = 

+5.85 


T 


(74) 


(75) 
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3.1.2 Calculation and Validation of Sensitivity Derivatives 


A study is conducted to validate the methodology which has been described herein for 
computing aerodynamic sensitivity derivatives, and to evaluate the efficiency of the procedures. 
In this study, geometric shape sensitivity derivatives are calculated using direct solution of Eq. 
(18) (together with some of the ancillary sensitivity equations presented in a previous sub- 
section), and are also computed using the method of “brute force” finite differences. In applying 
the method of finite differences, in order to insure great accuracy, central differences are used 
(requiring two CFD analyses per design variable) with an extremely small perturbation of each 
design variable, and the repeated CFD analyses are converged to machine zero. Sensitivity 
derivatives taken with respect to the ten previously described geometric shape design variables 
(Eqs. (65) through (75) ) are computed using these two methods, and are compared in Tables I, 
II, III, and IV. In Tables I and II. the sensitivity derivatives of C x and C y , respectively, for the 
lower wall of the nozzle are presented (C x and C y are aerodynamic force coefficients in the x 
and y directions, respectively, resulting from the integration of pressure and skin friction along 
the wall, and have been detailed previously.) In Tables III and IV, the sensitivity derivatives of 
the thrust, T. and the mass flow rate, m, of the nozzle are given. 


Design 

Variable 

X-Direction Force dC x 

Coefficient Sensitivity 

Direct 

Differentiation 

Finite 

Difference 

3 1 

-4.925 E+01 

-4.925 E+01 

$2 

-4.614 E+02 

-4.614 E+02 


+2.284 E+02 

+2.284 E+02 

3a 

-2.665 E+04 

-2.665E+04 

3s 

-8.327 E+01 

-8.327 E+01 

36 

-1.365 E-02 

-1.365 E-02 


+ 1.415 E+00 

+ 1.415 E+00 

3 $ 

P +6.235 E+00 

+6.235 E+00 

39 

-2.107 E+00 

-2.107 E+00 

3 io 

+3.886 E-01 

+3.886 E-01 


Table I - Comparison of the X-direction Force Coefficient, C x , Sensitivity 
Derivatives, Lower Wall, Double-Throat Nozzle, Laminar Flow 
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Design 

Variable 

Y-Direction Force 
Coefficient Sensitivity 

Direct 

Differentiation 

Finite 

Difference 

3i 

-3.024 E+02 

-3.024 E+02 

3 2 

+ 1.741 E+0 1 

+ 1.741 E+01 

h 

-2.625 E+01 

-2.625 E+01 


+ 1.664 E+03 

+ 1.664 E+03 

h 

+4.365 E-01 

+4.365 E-01 

■h 

+ 1.428 E+02 

+ 1.428 E+02 

■h 

-5.879 E+00 

^-5.879 E+00 

■h 

+2.331 E+02 

+2.331 E+02 

Jg 

-2.081 E+01 

-2.081 E+01 

^10 

+ 1.158 E+01 

+ 1.158 E+01 


Table II - Comparison of the Y-direction Force Coefficient, C y , 

Sensitivity Derivatives, Lower Wall, Double-Throat Nozzle, Laminar Flow 


Design 

Variable 

Thrust Sensitivity 

dT 

<Ok 

Direct 

Differentiation 

Finite 

Difference 

d, 

+ 1.841 E+02 

+ 1.841 E+02 


-4.869 E+00 

-4.869 E+00 

3 3 

+2.783 E+00 

+2.782 E+00 

3 4 

-3.254 E+02 

-3.254 E+02 

3s 

-7.304 E-01 

-7.303 E-01 

36 

^2.879 E-02 

-2.878 E-02 


-1.415 E+00 

-1.415 E+00 

3s 

-6.235 E+00 

-6.235 E+00 

3 9 

+2.106 E+00 

+2.106 E+00 

3 10 

-3.886 E-01 

-3.886 E-01 


Table III - Comparison of the Thrust Sensitivity 
Derivatives, Double-Throat Nozzle, Laminar Row 
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Design 

Variable 

Mass Flow Rate Sensitivity ^-jp- 



Direct 

Differentiation 

Finite 

Difference 

h 

+5.701 E+00 

+5.701 E+00 

■h 

-1.639 E-02 

-1.639 E-02 

h 

+2.819 E-02 

+2.819 E-02 


-2.168 E+00 

-2.168 E+00 

h 

-1.933 E-05 

r 1.895 E-05 

h 

-1.172 E-04 

-1.171 E-04 

h 

-1.096 E-09 

+0.000 E+00 


-4.158 E-09 

-5.684 E-09 


+ 1.131 E-04 

+ 1.134 E-04 

^10 

-8.740 E-13 

+0.000 E+00 


Table IV - Comparison of the Mass Flow Rate Sensitivity 
Derivatives, Double-Throat Nozzle, Laminar Row 


As expected, the agreement in the results above obtained using the two methods is excellent. 
Hie resulis only disagree for some of the mass flow rate sens, tiv, ties when these sensitivity 
denvauves are negligibly small. Table V is a comparison of die total CPU times which when, 
required to obtain one complete set of sensitivity derivatives. Clearly the direct differentiation 
approach is more efficient in this problem. However, due to die extreme care which was taken 
in insuring accuracy when using the finite difference approach, the CPU time reported here for 
this method should be taken as an upper bound “worst case”. 


CPU 

Direct 

Finite 

Time 

Differentiation 

Difference 


- 65 seconds 

- 1800 seconds 


Table V - Comparison of CPU times, Double-Throat 
Nozzle, Laminar Flow 


The direct differentiation approach of Eq. (18) has now been shown to efficiently and 
accurately yield aerodynamic sensitivity derivatives for 2D viscous internal flows which are 
laminar. Thts is the expected result, stnce the exact derivatives are found of the discrete algebraic 
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equations which model the fluid flow. Unfortunately, however, a major, ty of the practtcai flow 
problems of interest in aerodynamtc design are turbulent flows. If turbulence modeling is to be 
included in the flow calculations, then in principle, the equations which are used in the turbulence 
model which is selected should be consistently differentiated, and the results fully incorporated 
into the Jacobian matrices, [$] and [ff ] . of the fundamental sensitivity analysis / approximate 
analysis equations (i.e„ Eqs. (16) through (21)). In practice, of course, for typ.cal turbulence 
models, this represents a very difficult task, even impossible in the case of turbulence models 
which are not continuously differentiable (e.g.. the popular Baldwin-Lomax model). For this 
reason, in typical CFD codes, consistent linearization of the turbulence modeling is neglected 
in the construction of the implicit terms of the matrix, [|§], of Eq. (12) for the integral, on of 
the equations in time. Of course, the explicit spatial variation of Are turbulence modeling terms 
which is found in the residual vector on the right-hand side of Eq. (12) is also used on the 
lett-hand side in the terms ot' the Jacobian matrix, J£ 

As a potentially useful simplifying approximation in the calculation of sensitivity derivatives 
tor turbulent flows, an approach analogous to that described above for implicit time integration 
schemes is proposed. That is, differentiation of the turbulence modeling is neglected in the 
construction of the terms of the Jacobian matrices, [$], and [ff ] , of Eqs. (16) through (21). 
The explicit spatial variation of the steady-state values of the turbulence modeling terms is used, 
however, in both of these Jacobian matrices. This approximation is equivalent to a "locally 
constant” assumption applied to these terms. 

In order to judge the error in the sensitivity derivatives which might be expected in using this 
approximation, a preliminary investigation is conducted using the present double-throat nozzle 
problem, and the Baldwin-Lomax turbulence model [21j. Strictly speaking, this turbulence model 
is not applicable to the present low Reynolds number problem, but this does not conflict with 

the goals of this preliminary investigation into the effect of turbulence models on the accuracy 
of the sensitivity derivatives. 
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To begin the study, a converged steady-state numerical solution using the Baldwin-Lomax 
turbulence model is obtained for the previously described problem of laminar flow through the 
double-throat nozzle. Sensitivity derivatives are obtained for the turbulent case with respect 
to the same ten geometric shape design variables as for the laminar problem, again using the 
direct differentiation approach of Eq. (18), and also using the method of “brute force” finite 
differences. In using Eq. (18), the approximate treatment of the turbulence modeling terms is 
employed, as previously proposed. Tables VI and VII are a presentation and comparison of the 
computational results obtained using these two methods, where the sensitivity derivatives of the 
thrust and mass flow rate, respectively, are shown. 


Design 

Variable 

Thrust Sensitivity 

dT 

<Wk 

Direct 

Differentiation 

Finite 

Difference 

d, 

+ 1.823 E+02 

+ 1.821 E +02 

‘h 

-4.682 E+00 

-4.788 E+00 


+ 1.996 E+00 

+2.645 E+00 

3 4 

-2.961 E+02 

-3.140 E+02 

35 

-7.296 E-01 

-7.296 E-01 

3 6 

-5.707 E-02 

-2.470 E-02 

3y 

-1.402 E+00 

-1.475 E+00 


-3.074 E-00 

-3.624 E-00 

3 9 

+ 1.815 E +00 

+ 1.889 E+00 

3\o 

-2.627 E-01 

-2.705 E-01 


Table VI - Comparison of the Thrust Sensitivity 
Derivatives, Double-Throat Nozzle, Turbulent Flow 
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Design 

Variable 

Mass Flow Rate Sensitivity 

Direct 

Differentiation 

Finite 

Difference 

■h 

+5.703 E+00 

+5.697 E+00 


-1.699 E-02 

-1.641 E-02 

h 

+2.894 E-02 

+2.841 E-02 

h 

-2.238 E+00 

-2.187 E+00 

■is 

-2.317 E-05 

-2.463 E-05 


-6.753 E-04 

-1.516 E-04 

■h 

-5.658 E-10 

+0.000 E+00 

'h 

-1.949 E-09 

-4.263 E-09 

^9 

+2.816 E-04 

+ 1.045 E-04 

•-?10 

-2.084 E-13 

-1.421 E-09 


Table VII - Comparison of the Mass Flow Rate Sensitivity 
Derivatives, Double-Throat Nozzle, Turbulent Flow 


Note that in generating the brute force” sensitivity derivatives shown in Table VI and 
VII, the same extreme care was taken to insure great accuracy of the terms that was taken 
previously for the laminar results. Therefore, the “brute force" values shown for the turbulent 
case are considered to be more accurate than the values obtained using the method of direct 
dilterentiation (because of the approximate treatment of the turbulence modeling terms which 
was used in the latter approach). In general, however, the agreement between the two methods 
is reasonably good, hopefully good enough for future use in a design optimization strategy. 

3.1.3 Optimization Problem — Double-Throat Nozzle 

Having verified the accuracy and efficiency of the methodology for computing sensitivity 
denvauves, a short study is conducted to demonstrate the application of the sensitivity derivatives 
in a model optimization example problem, using the double-throat nozzle described above. A 
generalized procedure is proposed for aerodynamic design optimization, where the principal 
computational tasks are sub-divided into distinct modules, or major subroutines, all controlled 
by a central, external, general purpose design optimization computer program. The well-known 
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design optimization program selected for use in the present study is called ADS (Automated 
Design Synthesis), which is documented in Ref. [48], although others are available. 

The principal computational tasks to be performed dynamically in the automated aerodynamic 
design strategy are accomplished using a standard CFD code (in order to re-evaluate the 
aerodynamic performance variables), a sensitivity analysis program, based on the CFD analysis 
code (lor computing the sensitivity derivatives required by the optimization program), a grid 
regeneration program (for shape optimization, as the shape of the domain changes) and grid 
sensitivity capability (necessary in computing shape sensitivity derivatives). In addition, an 
optional approximate analysis capability (based on Eqs. (19) or (21), for geometric shape 
changes) can be included as a subset of the sensitivity analysis module, in order to provide 
significantly more efficient re-evaluation of the aerodynamic performance variables, when the 
changes in the design variables are small (i.e., less than some specified tolerance). The validity, 
accuracy, and efficiency of the approximate analysis methodology is documented in Refs. 
[3,5,6,7]. The proposed aerodynamic design optimization strategy is illustrated in Fig. (10). 

Starting with the initial nozzle shape and steady-state solution, thv. proposed aerodynamic 
shape optimization method of Fig. (10) is applied to the nozzle problem, as follows: 

1) The force coefficient in the x direction, C x , (henceforth known as the objective function), 
is to be minimized in the design, subject to the explicit constraints on the design variables, 
given subsequently. C x represents the internal aerodynamic drag on the lower wall of 
the nozzle. 

2) The design variables are the ten elements of 3 , defined in Eq. (75), and an upper and 
lower bound is placed on each, where the total maximum allowable change in each 
during the design, either an increase or decrease, is specified a priori to be no greater 
than 15% of the initial value of each variable. 

3) Sensitivity derivatives with respect to the design variables of the objective function are 
obtained using the method of direct differentiation (Eq. (18)) and the ancillary sensitivity 
equations given previously. Sensitivity derivatives of the explicit constraints (i.e., the 
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upper and lower bounds) on the design variables are not needed, as they do not change 
during optimization. 

4) Mesh regeneration and grid sensitivity derivatives are provided by the grid generation 
computer program wh.ch produced the initial mesh, and gnd sensittvity is obtained using 
this piogram and the “brute force” approach (i.e., Eq. (37)). 

Following application of the automated destgn optimtzalton strategy, the results are sum- 
marized in Tables VIII and IX. 


Performance 

Variable 

Initial 

Design 

Final 

Design 

Percent 

Change 

Force Coefficient. C x 
(Objective Function) 

+754.4 

+313.9 

-58.39 

Thrust, T 

+ 190.4 

+205.4 

+7.88 

Mass Flow Rate, rii 

+5.463 

+6.284 

+15.0 

Specific Thrust, F s 

+34.85 

+32.69 

-6.20 


Table VIII Comparison of the Aerodynamic Performance 
of the Initial and Final Designs, Nozzle Problem 


Design 

Initial 

Final 

Percent 

Variable 

Design 

Design 

Change 

ijjL 

+ 1.0 

+ 1.150 

+15.00 

ih 

-4.0 

-3.400 

+15.00 

3s 

+2.3 

+2.017 

-12.30 

3 4 

-0.03 

^ -0.0255 

+15.00 

3s 

-10.0 

-8.500 

1" + 15.00 

3 6 

-12.0 

12.000 

0.0 

$7 

+7.0 

+7.000 

0.0 

^8 

+ 1.6 

+ 1.599 

-0.0625 

3 Q 

+ 14.0 

+ 14.000 

0.0 

3 io 

+5.85 

5.850 

0.0 ' 


Table IX - Comparison of the Initial and Final Design Variables. Nozzle Problem 
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Clearly a significant reduction in the aerodynamic drag on the internal walls of the nozzle (as 
defined in the present problem statement) was achieved, as demonstrated by the desired reduction 
in the objective function in the new design. Figs. (11a) and (lib) represent a comparison of 
the geometric shapes of the initial and final nozzle designs, respectively. (Fig. (11a) is a repeat 
of Fig. (7).) In reducing the aerodynamic drag on the internal walls of the nozzle, the shape 
of the final design clearly displays a predictable pattern of “smoothing out the bumps” in the 
contours of the walls. Had the explicit bounds on the design variables been removed, this 
effect would surely have been even greater. The results shown here were obtained following 
only three evaluations of the objective function (following the initial evaluation). The option of 
approximate analysis was not used in the present test problem. 

In the present model design problem, aerodynamic drag reduction is considered, as it is 
of general interest in many aerodynamic design problems. In the design of nozzles, however, 
clearly the thrust, mass flow rate, and specific thrust, F s , (F s = T/m) are variables of special 
interest to the designer. For this reason, the values of these variables for the initial and final 
designs have been included in Table VIII. A more realistic shape optimization problem for 
nozzles might typically select the maximization of F s as the objective function, perhaps with 
the side constraint that the thrust of the final design be no less than that of the initial design. 
Reformulation of the present test problem in this manner is straightforward, in principle, since 
the sensitivity derivatives of these terms are available (i.e„ thrust and mass flow rate sensitivities 
have been presented herein, from which specific thrust sensitivities are obtained using Eq. (64). 

3.2 External Flow — NACA 4-Digit Airfoil Problem 
3.2.1 Description of the Test Problem 

The external flow problem which is studied in the present work is that of an isolated airfoil. 
The initial airfoil which is selected is the NACA 2412, the profile of which is defined by Eqs. 
(41) and (42), with maximum thickness, T = 0.12, maximum camber, C = 0.02, and location of 
maximum camber, L = 0.40. The computations are performed on a “C” mesh using 128 points 
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in the “around the airfoil” direction, (with 68 of these on the surface of the airfoil) and 32 points 
in the normal to the airfoil direction, with grid stretching near the airfoil surface. The tar-held 
boundary is placed twenty chords from the airfoil. The grid is generated using the algebraic mesh 
generation program of Ref. [42], The gnd used in the present test case is clearly inadequate 
for use in resolving the full physics of viscous flow over an airfoil. Therefore, it should be 
understood that the objective of the work for which results are to be presented is not at this point 
to produce highly accurate calculations of the flow held. The objective rather is to successfully 
implement, demonstrate, and validate the ideas on a computationally inexpensive model problem. 

A converged conventional numerical solution is obtained for the initial geometry and grid for 
laminar flow, for a freestream Mach number, M*, = 0.70, Reynolds number, RE L = 5000.. and 
angle of attack a = 0.0°. Boundary conditions on the airfoil surface, across the “wake-cut”, and 
at the tar-held boundaries have been discussed previously. When the “lift corrected” far-held 
boundary conditions are not used, the lift <C L >, drag (C D ), and moment (C M ) coefficients for the 
calculation were C L = 0.1232, C D = 0.06824, and C M = - 0.05328. When the lift correction 
of Ret. [40] js applied at the tar-held boundaries, these coefficients become Q = 0.1252, Co 
= 0.06821, and C M = - 0.05325. Therefore, the numerical effect of the far-held lift correction 
is not strong in this example problem, which is attributed to the rather large extent of the far- 
held boundary from the airfoil (twenty chords) which was used. This effect will become more 
significant, however, as the far-held boundary is brought closer to the airfoil surface [40], i.e.. 
the lift-corrected boundary conditions become more necessary. Of course, the freedom to move 
the far-held boundary closer to the airfoil surface allows more effective use of the grid points 
which are available, in order to more accurately resolve the flow physics. 

3.2.2 Calculation and Validation of Sensitivity Derivatives 

A study is conducted to validate the methodology which has been described herein for 
computing aerodynamic sensitivity derivatives for airfoils, and to evaluate the efficiency of 
the procedures. Unless otherwise stated, sensitivity derivatives are computed using the steady- 
state solution without the point-vortex lift-corrected far-held boundary conditions. In particular. 
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geometric shape sensitivity derivatives are calculated by solution of Eq. (18) using the hybrid 
direct solver / iterative strategy proposed in Eq. (30) (and also using some of the ancillary 
sensitivity equations presented previously). For comparison purposes, the sensitivity derivatives 
are also computed using the ‘brute torce ’ method of finite differences. In applying this approach, 
to insure good accuracy, central finite difference approximations are used with a forward and 
backward perturbation of each design variable of 0.01%, and the repeated CFD analyses are 
converged to machine zero. 

The sensitivity derivatives ot the lilt, drag, and moment coefficients are computed with 
respect to the three geometric shape parameters, T, C, and L of Eqs. (41) and (42) using 
the initial airfoil, grid, and flow conditions which have been described immediately above. 
In solving Eq. (18), the required grid sensitivity terms (i.e„ j||-}) were provided using 
the previously described “elastic membrane” representation of the computational domain, and 
the prescribed displacement approach, together with Eqs. (39) through (45). In addition, 
as a consistency requirement, this grid methodology is also used in generating the necessary 
perturbed grids which are required in computing aerodynamic sensitivity derivatives by the 
“brute force” approach. 

When solving Eq. (18) using the hybrid direct solver / iterative method of Eq. (30), it 
is necessary to establish some type of criterion by which the iterative solution of the linear 
system tor each design variable is judged to be “converged,” and is terminated. To study this, 
the sensitivity derivatives are calculated with Eq. (30) using progressively stricter convergence 
criteria, starting with a calculation where no iterations are used (i.e., the “periodic” terms outside 
of the central bandwidth are simply neglected) and ending with a calculation where the average 
global error is reduced four orders-of-magnitude. A summary of these calculations, including 
the results of the finite difference approach, is presented in Tables X, XI. and XD, one table for 
each of the three geometric shape design variables, T, C, and L, respectively. 

As expected, as the reduction in the error by the iterative process defined by Eq. (30) is 
increased, the agreement between the sensitivity derivatives calculated using the method of direct 
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differentiation (Eq. (18)) and the same terms calculated using the method of finite differences 
is excellent. In addition, the validity of the “elastic membrane” methodology for computing 
grid sensitivities is confirmed. Furthermore, in the present problem, it is demonstrated that 1) 
approximately a three orders-of-magnitude reduction in the error is more than satisfactory for 
engineering purposes in applying Eq. (30), but 2) complete neglect of the implicit periodic 
boundary condition terms is completely unacceptable in using Eq. (18) to predict the sensitivity 

derivatives (Note: In contrast, these terms typically are neglected in the integration of the 
equations in time, using Eq. (12)). 


Method 

of 

Solution 

Lift Sensitivity 

dOu - d£, 

d U ~ ~3T 

Drag Sensitivity 

XV, XV 

~ dT 

Moment Sensitivity 

3Cm dCu 

diii ~ dT 

No. of 
Iterations 
Using Eq. (30) 

Eq. (30), 

Without 

Iterations 

-9.334 E-01 

+4.723 E-01 

+9.815 E-02 

1 

Eq. (30), 

1 OM 

-2.859 E+00 

+4.267 E-01 

+4.755 E-01 

13 

Eq. (30), 
2 OM 

-3.117 E+00 

+3.972 E-01 

I — 

+5.278 E-01 

64 

Eq. (30), 
3 OM 

-3.126 E+00 

+3.939 E-01 

+5.307 E-01 

219 

Eq. (30), 

4 OM 

-3.126 E+00 

+3.938 E-01 

+5.307 E-01 

300 

Finite 

Difference 

-3.126 E+00 

+3.938 E-01 

+5.307 E-01 

N/A 


*OM Reters to the number of Orders-of-Magnitude reduction in the average global error. 

Table X - Comparison ot Lift and Drag Sensitivity Derivatives 
With Respect To Maximum Thickness, tfi=T, NACA 2412 Airfoil 
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Method 

of 

Solution 

Lift Sensitivity 

XV ,xy 

~ x 

Drag Sensitivity 

Xn X'n 

SJj ~ SC 

Moment Sensitivity 

ssf ~ -sc 1 

No. of 
Iterations 
Using Eq. (30) 

Eq. (30), 

Without 

Iterations 

+5.206 E+00 

+3.429 E-01 

-2.520 E+00 

1 

Eq. (30), 

l OM 

+4.175 E+00 

+3.780 E-01 

-2.261 E+00 

14 

Eq. (30), 
2 OM 

+3.988 E+00 

+3.663 E-01 

-2.225 E+00 

39 

Eq. (30), 
3 OM 

-3.968 E+00 

+3.603 E-01 

-2.220 E+00 

188 

Eq. (30), 
4 OM 

+3.968 E+00 

+3.603 E-01 

-2.220 E+00 

^76 

Finite 

Difference 

+3.968 E+00 

+3.603 E-01 

-2.220 E+00 

N/A 


*0M Reters t0 the number of Orders-of-Magnitude reduction in the average global error. 

Table XI - Comparison of Lift and Drag Sensitivity Derivatives 
With Respect To Maximum Camber, ^ 2 =C, NACA 2412 Airfoil 


Method 

of 

Solution 

Lift Sensitivity 

x ,, x, 

X, SL 

Drag Sensitivity 

Xn _ X^ 

X 3 ~ SL 

Moment Sensitivity 

x M _ ac M 

X 3 SL 

No. of Iterations 1 
Using Eq. (30) 

Eq. (30), 

Without 

Iterations 

-4.293 E-02 

-3.899 E-03 

-3.865 E-02 

1 

Eq. (30), 
1 OM 

-1.466 E-02 

-3.422 E-03 

-4.456 E-02 

5 

Eq. (30), 
2 OM 

-1.869 E-02 

-3.334 E-03 

-4.386 E-02 

24 

Eq. (30), 
3 OM 

-1.819 E-02 

-3.304 E-03 

-4.396 E-02 

49 

Eq. (30), 
4 OM 

-1.816 E-02 

-3.290 E-03 

-4.398 E-02 

195 

Finite 

Difference 

-1.815 E-02 

-3.290 E-03 

-4.398 E-02 

N/A 


*OM Refers to the number of Orders-of-Magnitude reduction in the average global error. 

Table XII - Comparison of Lift and Drag Sensitivity Derivadves With 
Respect To Location of Maximum Camber, /i 3 =L, NACA 2412 Airfoil 
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Table XIII is a comparison of the relative computational cost of obtaining the sensitivity 
derivatives using the procedures for which results have been given in the three preceding tables. 
In the present problem, clearly the finite difference method is less efficient than the use of Eq. 
(30). even when the required error reductions using Eq. (30) are specified to be large. However, 
because ot the great care which was taken in applying the finite difference method to insure 
very accurate results (for validation purposes), the CPU times reported for this finite difference 
method should be taken to represent a "worst case" scenario. 


Method of 
Solution 

CPU Time, Total 
(seconds) 

Eq. (30), Without 
Iterations 

27 

Eq. (30) 
1 OM 

33 

Eq. (30) 
2 OM 

54 

Eq. (30) 
3 OM 

124 

Eq. (30) 

4 OM 

191 

Finite 

Difference 

840 


*OM Refers to the number of Orders-of-Magnitude reduction in the average global error. 

Table XIII - Comparison of CPU Times - Airfoil Problem 

As mentioned earlier, in applying the hybrid direct solver / iterative strategy of Eq. (30) 
in the solution ot Eq. (18), for the results presented above, it was necessary to use under- 
relaxation to force the iterations to converge. A single under-relaxation parameter, w = 0.75 
was used. The required use of this parameter, in addition to the relatively large numbers of 
iterations (and hence the proportionally large amounts of CPU time) required in reducing the 
error to an acceptable level is worrisome, and was neither the expected nor the desired result. 
As stated earlier, these difficulties were not encountered in the related work of Ref. [33,39], 
Additional work is needed to correct these difficulties, and to improve the efficiency of the 
algorithms for solving the equations of sensitivity analysis. Recent studies have cast these large 
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systems of linear equations in "delta” or “incremental” form [49,50], where it is seem that these 
difficulties with matrix ill-conditioning are overcome. This incremental formulation allows for 
the accurate efficient iterative solution of these equations using the identical approximate left- 
hand side coefficient matrix operator and algorithms which are commonly using in the numerical 
time integration of the non-linear flow equations to steady-state [26,27,34,35,36,37.38], 

For airfoil design problems, it the sensitivity derivatives of only two or three performance 
variables (i.e„ C L , C D , and maybe Cm) are needed, but the number of design variables is 
significantly larger than two or three, then the issue of a slow convergence rates during multiple 
iterative solutions of Eq. (30) can be addressed to a large extent by switching from the present 
method to the adjoint-variable method for computing the sensitivity derivatives [2,4], The 
computational work of the adjoint-vanable approach is completely independent of the number of 
design variables, and requires the solution (either directly or iteratively) of a single large system 
ot linear equations (having coefficient matrix, -[§§] T ) for each aerodynamic performance 
variable tor which the sensitivity derivatives are required, which as stated is typically only two 
or three variables, for airfoils. Clearly this is potentially a very large savings if the convergence 
rate ot Eq. (30) is slow and the number of design variables is large. 

The sensitivity derivatives of the steady-state solution with the lift-corrected far-field bound- 
ary conditions [40] was considered next, using the methodology proposed in Eqs. (31) through 
(33). The mo dified hybrid direct solver/conventional iterative solver strategy of Eqs. (35) (in- 
cluding Eq. (36)) was applied in the solution for the sensitivity derivatives. The results are 
summarized in Table XIV, where it is seen that the agreement in the calculations is very good in 
comparing the results obtained using Eqs. (35) and (36) (following a four orders-of-magnitude 
reduction in the average global error) with the results obtained using the method of “brute 
torce finite differences. Of additional interest is to compare the results of Table XIV with the 
corresponding results presented in Tables X, XI, and XII (i.e., the comparison is of the sensi- 
tivity derivatives of the solution with the lift-corrected boundary conditions with the sensitivity 
derivatives ot the solution without this correction). 
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Solution 

Method 

Design 

Variable 

dC' L 

dC D 

dC\i 

dT.C.L 

dT. C. L 

dT.C.L 

Eqs (35) and (36), 
4 OM* 

T 

-3.211 E+00 

+3.959 E-01 

+5.337 E-01 

C 

+4.015 E+00 

+3.578 E-01 

-2.216 E+00 

L 

-1.751 E-02 

-3.327 E-03 

-4.410 E-02 

"Brute Force." 
Finite Difference 
Method. 

T 

-3.211 E+00 

+3.959 E-01 

+5.337 E-01 

C 

+4.015 E+00 

+3.578 E-01 

-2.217 E+00 

L 

-1.752 E-02 

-3.326 E-03 

-4.410 E-02 


*0M Refers to the number ot Orders-of-Magnitude reduction in the average global error 
Table XIV — Comparison ot the Sensitivity Derivatives With The 
Lift-Corrected Far-Field Boundary Conditions - Airfoil Problem 


3.2.3 Optimization Problem — NACA 4-digit Airfoil 

Having successfully validated the accuracy and efficiency of the proposed strategies for 
computing grid regeneration, grid sensitivity derivatives, and finally, aerodynamic sensitivity 
derivatives tor isolated airfoils on “C” and “O” meshes, the application of these methods is 
demonstrated here in a model airfoil design optimization problem. The optimization strategy 
described previously tor the nozzle test problem (and illustrated in Fig. (10)) is applied, using 
the NACA 2412 airtoil, grid, and steady-state solution described above as the initial design. The 
optimization problem is formulated as follows: 

1) The lift coefficient, Cl, is to be maximized in response to geometric shape variations, 
subject to the design constraints defined subsequently. 

2) The drag coefficient, Cd. is to be no greater than the value of the initial design. 

3) The design variables (i.e., the elements of /3) are the three parameters of the NACA 
4-digit airfoil, i.e., T, C, and L, defined previously. The “linearized” NACA 4-digit 
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airtoil ot Eqs. (43) and (44) defines the range of allowable airfoil shapes in terms of 
these three parameters. (As stated previously, this linearized equation only approaches 
the true NACA 4-digit parameterization of Eqs. (41) and (42) for small changes in 
the parameters, but the current parameterization is not limited to small changes). In 

addition, upper and lower bounds are defined a priori, for each of the design variables, 
as follows: 

0.08 < T < 0.16 
0.00 < C < 0.05 
0.20 < L < 0.60 

4) Sensitivity derivatives of the objective function, C L , and the constraint. C D , with respect 
to the design variables are obtain by solution of Eq. (18), using the hybrid direct solver 
/ iterative strategy of Eq. (30). The global average error reduction in the solution of Eq. 
(18) is specified to be three orders-of-magnitude for each design variable. Sensitivity 
derivatives of the explicit constraints on the design variables (i.e., the upper and lower 
bounds) are not needed. 

5) Mesh regeneration and grid sensitivity derivatives are provided dynamically during the 
automated optimization process using the finite element method and the “elastic mem- 
brane” prescribed displacement model of the computational grid, presented previously. 

One important benefit and use of sensitivity analysis is the ability to study the sensitivity 
denvauves of the initial design before the design process begins, in order to identify which 
design variables are important (i.e., very sensitive) to the problem, and which are not. By 

this, often it can be safely judged that some of the design variables are unneeded, and can be 
eliminated from the problem. 

As an example of this, examination of the sensitivity derivatives in Tables X, XI, and XII (or 
Table XIV) reveals that the lift and drag sensitivities are about two orders-of-magnitude smaller 
for the third design vanable, L, compared to the sensitivity derivatives for the remaining two. 


50 



This means that the third design variable is not important to the current design problem, and 
for increased efficiency in the design, can be eliminated. To verify this, the design optimization 
methodology was implemented twice, once including the design variable. L. and once neglecting 
it completely (i.e., fixing it to be the value of the initial design). As expected, there was no 
significant difference following optimization in the final improved designs, either in aerodynamic 
performance, or in the final values ol the remaining two design variables, when the third 
design variable was discarded from the problem. The results of the airfoil optimization shown 
subsequendy in Tables XV and XVI are therefore for the case where the third design variable. 
L, has been dropped from the design problem. 


Performance 

Initial 

Improved 

Percent 

Variable 

Design 

Design 

Change 

Cl, Lift Coefficeint 

0.1232 

0.3658 

+ 196.9 

(objective function) 




Cd, Drag Coefficient 

0.06824 

0.06791 

-0.48 

(a constraint) 




Lift/Drag Ratio 


5.387 

+ 198.4 


Table XV - Comparison of the Aerodynamic Performance 
of the Initial and Final Designs, Airfoil Problem 


Design 

Initial 

Improved 

Percent 

Variable 

Design 

Design 

Change 

tfi=T 

0.1200 

0.0800 

-33.33 

3 2 =C 

0.0200 

0.0455 

+127.5 

3 }= L 

0.4000 

0.4000 

0.0 


Table XVI - Comparison of the Initial and Final Designs Variables, Airfoil Problem 

Clearly an airfoil shape having significantly improved aerodynamic performance for the 
given conditions has resulted from the automated design methodology, in the present test case, 
as demonstrated by the increased lift coefficient and lift / drag ratio of the new design. Figs. (12a) 
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and (12b) represent a comparison between the geometric shapes of the initial and final airfoil 
designs, respectively. The results shown here were obtained following only two evaluations of 
the objective function (following the initial evaluation). The option of approximate analysis was 
not used in the present test problem. 

4.0 Summary and Conclusions 

A gradient based design optimization strategy has been developed in detail for use in practical 
aerodynamic design problems, using the 2D thin-layer Navier-Stokes equations. The strategy 
is general in nature, and is based on the classic idea of constructing different modules for 
performing the major tasks such as function evaluation, function approximation and sensitivity 
analysis, mesh regeneration and grid sensitivity analysis, all driven and controlled by a general 
purpose design optimization program. The methodology has been successfully demonstrated on 
both an internal and an external flow problem, where in each case, a significant improvement 
in aerodynamic performance was achieved. 

In developing the methods which have been presented herein, the major effort has been fo- 
cused on techniques tor efficiently and accurately calculating the necessary sensitivity derivatives 
tor use by the optimization program. In particular, much effort was spent in successfully devel- 
oping and implementing a consistent treatment of the boundary conditions in the calculation of 
the aerodynamic sensitivity derivatives for the classic problem of external flow over an isolated 
lifting airfoil on “C” or “O” meshes. As a minor extension of the methodology for computing 
the sensitivity derivatives, an efficient strategy for approximate analysis is available, a capability 
which is potentially very useful in a design environment. 

One of the most difficult design variables to consider is that of geometric shape, which 
has been the emphasis of the present work. With geometric shape, the issue of grid sensitivity 
analysis and grid regeneration must be considered in the development of an automated gradient 
based optimization scheme. In the present work, an efficient method for handling this difficulty 
based on well-known ideas from the discipline of structural shape design optimization has been 
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successfully applied to the aerodynamic design problem. The technique employs the hnite 
element method, and an elastic membrane representation of the computational domain. 

The investigation which has been presented herein is still in its initial stages, and there is 
an unlimited amount left to be done. For example, for the viscous airfoil problem which was 
presented herein, the grid was of insufficient size to truly resolve the full physics of a viscous 
flow. Future work will focus on implementation of the methods demonstrated here to full scale 
grids. A more complete parameterization of the airfoil boundary with more design variables is to 
be used, in order to increase the range of the design space. Design with turbulence modeling will 
be attempted. Finally, development of more efficient strategies (with respect to computational 
work and computer storage) for computing the sensitivity derivatives is underway. 
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A Symbol Indicates Simple Support, Zero Displacement. 

(Points Supported Included Leading And Trailing Edges, 
Points In Wake, And All Points At The Far-Field Boundary) 


Fig. (3) - Illustration Of The “Elastic Membrane” Representation Of The Computational 
Domain, Using The Fictitious Load Method For Computing Grid Sensitivity Derivatives. 
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A Symbol Indicates Simple Support, Zero Displacement. 
(Points Supported Included Leading And Trailing Edges, 
Points In Wake, And All Points At The Far-Field Boundary) 


Fig. (4) - Illustration Of The “Elastic Membrane” Representation 
Of The Computational Domain Using The Prescribed Boundary 
Displacement Method For Computing Grid Sensitivity Derivatives. 
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Fig. (9) Illustration Of The Parameterization Of The 
Boundary Shape, Double-Throat Nozzle Problem. 



Fig. (IQ - Illustration of The Flow Chan For The Proposed 
Sensitivity Based Aerodynamic Shape Optimization Strategy. 









Fig. (12a) - Illustration Of The Initial Airfoil Design, 
NACA 2412, Moo = 0.70, RE L =5000, a=0.0° 


Fig. ( 12b) - Illustration Of The Final Airfoil Design 
Following Optimization, Moo =0.70, RE L =5000, a=0.0° 


